#include "mecafluide.hpp"

#include <research/life-workspace/Fluid-Structure/functionSup.cpp>

//#include <mecafluide_cl.cpp>

namespace Feel
{
namespace FSI
{

    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    void
    NavierStokes<Dim, Order, Entity>::updateResidual( const vector_ptrtype& X, vector_ptrtype& R )
    {

        using namespace Feel::vf;

        if (M_LookComment) std::cout <<"[NavierStokes] : updateResidual start \n";

        boost::timer time;
        time.restart();

        this->changeRepository();

        //--------------------------------------------------------------------------------------------------//

        auto mesh = this->mesh();
        auto Xh = this->functionSpace();

        element_fluid_type U( Xh, "u" ); U = *X;
        element_fluid_type V( Xh, "v" );
        element_fluid_0_type u = U.template element<0>();
        element_fluid_0_type v = V.template element<0>();
        element_fluid_1_type p = U.template element<1>();
        element_fluid_1_type q = V.template element<1>();


        //--------------------------------------------------------------------------------------------------//

        //Tenseur des deformations
        auto defv = 0.5*(gradv(u)+trans(gradv(u)) );
        auto def = 0.5*(grad(v)+trans(grad(v)) );
        //Matrice Identitee
        auto Id = oneX()*trans(oneX()) + oneY()*trans(oneY());
        // Tenseur des contraintes (trial)
        auto Sigmav = -idv(p)*Id + 2*M_mu*defv;
        // Tenseur des contraintes (test)
        auto Sigma = -id(q)*Id + 2*M_mu*def;
        //force volumique :
        auto f = MECAFLUIDE_FORCE_VOLUMIQUE(M_mu);

        auto clDef = MECAFLUIDE_CL(M_time);
        //--------------------------------------------------------------------------------------------------//


        double EE = 3*1.0e6;//21*1e5;
        double TTsigma = 0.3;//45;//0.3;//0.28;
        //EE = 0.2;//module young (caoutchou)
        //sigma = 0.5;//coeff de poisson (caoutchou)
        double TTmu = EE/(2*(1+TTsigma));
        double TTlambda = EE*TTsigma/((1+TTsigma)*(1-2*TTsigma));

        element_fluid_velocity_type displStrucRef((M_Xh_ref->template functionSpace<0>()));
        displStrucRef = vf::project(displStrucRef.functionSpace(),elements(displStrucRef.functionSpace()->mesh()),idv(*M_aux));
        element_fluid_velocity_type displStruc((Xh->template functionSpace<0>()));

        //std::cout << "\n Nb dof : " << displStruc.functionSpace()->nDof();
        //std::cout << "\n Nb dof : " << M_aux->functionSpace()->nDof();

        displStruc.container()=displStrucRef.container();
        auto Structepsv = 0.5*(gradv(displStruc)+trans(gradv(displStruc)));
        auto Structsigmav = TTlambda*trace(Structepsv)*Id + 2*TTmu*Structepsv;

        //--------------------------------------------------------------------------------------------------//

        form1( Xh, R, _init=true );

        bool M_isNLTerme=true;
        if (M_isNLTerme)
            if (M_isMoveDomain)
                {
                    //le terme non linéaire
                    form1( Xh, R ) +=
                        integrate( elements(mesh), _Q<3*(Order)>(),
                                   //M_rho*trans(idv(u))*gradv(u)*id(v) );
                                   M_rho*trans(idv(u) -idv(*M_mesh_velocity))*(gradv(u))*id(v) );
                                   //M_rho*trans(id(v))*(gradv(u))*(idv(u) - idv(*M_mesh_velocity)) );
                }
            else
                {
                    //le terme non linéaire
                    form1( Xh, R ) +=
                        integrate( elements(mesh), _Q<3*(Order)>(),
                                   M_rho*trans(id(v))*(gradv(u))*(idv(u) ) );
                }
#if 0
        // sigma : grad(v) sur Omega
        form1( Xh, R ) +=
            integrate( elements(mesh), _Q<2*(Order)>(),
                       + trace(Sigmav*trans(grad(v))) );
#else
        form1( Xh, R ) +=
            integrate( elements(mesh),
                       + 2*M_mu*trace(trans(defv)*grad(v)) );
                       //+ 2*M_mu*trace(defv*trans(grad(v))) );
                       //+ 2*M_mu*trace(trans(defv)*def) );
        form1( Xh, R ) +=
            integrate( elements(mesh),
                       -idv(p)*div(v) );
        /*
        form1( Xh, R ) +=
            integrate( boundaryfaces(mesh),
                       - idv(p)*trans(id(v))*N() );
        */
#endif

#if 1
        form1( Xh, R ) +=
            integrate( elements(mesh),// _Q<2*(Order)>(),
                       - M_rho*divv(*M_mesh_velocity)*trans(idv(u))*id(v) );
#endif


        // terme assurant l'incompressibilité
        form1( Xh, R ) +=
            integrate( elements(mesh), _Q<2*(Order-1)>(),
                       + divv(u)*id(q) );


        ForEachCL( clDef,cl::neumann_scal,
                   form1( Xh, R ) +=
                   integrate( markedfaces(mesh,mesh->markerName(PhysicalName)),
                              - trans (Expression*N())*id(v) )     );

        //--------------------------------------------------------------------------------------------------//

        if (M_weakCL)
            {

                ForEachCL( clDef, cl::dirichlet_vec,
                           form1( Xh, R ) +=
                           integrate( markedfaces(mesh,mesh->markerName(PhysicalName)),// _Q<2*Order>(),
                                      //- trans(Sigmav*N())*id(v)
                                      - 2*M_mu*trans(defv*N())*id(v) + trans(idv(p)*N())*id(v)
                                      + penalbc*trans(idv(u)/*-Expression*/)*id(v)/hFace() ) );

                ForEachCL( clDef, cl::paroi_mobile,
                           form1( Xh, R ) +=
                           integrate( markedfaces(mesh,mesh->markerName(PhysicalName)),// _Q<2*Order>(),
                                      //- trans(Structsigmav*N())*id(v)
                                      - trans(Sigmav*N())*id(v)
                                      //- trans(Sigma*N())*idv(u) // sym
                                      //- 2*M_mu*trans(defv*N())*id(v) + trans(idv(p)*N())*id(v)
                                      + penalbc*trans(idv(u)/*-idv(this->meshVelocity())*/)*id(v)/hFace() ) );

            }
        //--------------------------------------------------------------------------------------------------//

        //Force Volumique
        form1( Xh, R ) +=
            integrate( elements(mesh), _Q<Order>(),
                       - trans(f)*id(v) );




        //------------------------------------------------------------------------------------//

        if (!M_isStat)
            {

                form1( Xh, R ) +=
                    integrate( elements(mesh), _Q< 2*Order >(),
                               + M_rho*trans(idv(u))*id(v)*M_bdf_fluid->polyDerivCoefficient(0) );

                auto Buzz = M_bdf_fluid->polyDeriv();
                auto buzz = Buzz.template element<0>();

                form1( Xh, R ) +=
                    integrate( elements(mesh), _Q< 2*Order >(),
                               - M_rho*trans(idv(buzz))*id(v) );
            }

        //------------------------------------------------------------------------------------//

        double penalisationBeta=0.05;
        auto beta_N = val(trans(N())*M_rho*(idv(u)-idv(*M_mesh_velocity)));
        double order_scaling = math::pow( double(Order), 3.5 );
        form1( Xh, R ) +=
            integrate( internalfaces(Xh->mesh()), //typename MyIm< 4*uOrder-2 >::type(),
                       val(penalisationBeta + vf::pow(hFace(),2.0)*( beta_N*beta_N ))*(trans(jumpv(gradv(u)))*jump(grad(v)) )/order_scaling );

#if 0
        int stabilisation = this->application()->vm()["stab"].template as<int>();
        if (stabilisation==1) {
            //terme de stabilisation
            form1( Xh, R ) +=
                integrate( internalfaces(mesh),// _Q< 3*Order >(),
                           // + 0.025*vf::pow(hFace(),2.0)*( abs( trans(idv(u))*N() ))*(trans(jumpv(gradv(u)))*jump(grad(v)) )  );
                           + (0.05+vf::pow(hFace(),2.0))*( (trans(idv(u))*N())*( trans(idv(u))*N() ))*(trans(jumpv(gradv(u)))*jump(grad(v)) )  );
        }
#endif

        //------------------------------------------------------------------------------------//

        R->close();

        if (true)//!M_weakCL)
            {
                auto dirElem = this->getDirichletElement();
                ForEachCL( clDef,cl::dirichlet_vec,
                           modifVec222(*dirElem,R, vf::vec(vf::cst(0.),vf::cst(0.)) , PhysicalName ) );
                ForEachCL( clDef,cl::paroi_mobile,
                           modifVec222(*dirElem,R, vf::vec(vf::cst(0.),vf::cst(0.)) , PhysicalName ) );
            }


        double ti = time.elapsed();
        if (M_LookComment) std::cout <<"[NavierStokes] : updateResidual finish in "<<ti<<"s\n";



    }

} // end namespace NavierStokes
} // end namespaceFeel


// instantiation
template class Feel::FSI::NavierStokes<MECAFLUIDE_DIM,MECAFLUIDE_ORDER_VELOCITY,Feel::Simplex>;
//template class Feel::NavierStokes<2,4,Feel::Simplex>;
